
% this is the analog of poolStiMData_v1_2 but for pooling data from diving
% vessels;
%
% inputs:
% regexp is the regular expression that preceeds the movie # for all the
% movies in the directory
%
% outputs:
% trajectoriesDDD - each column is a time point aligned to the start of the
% stimulus (offset by prestimTime, fps, poststimTime) in deltaD/D array;
% baseline is determined as average position over the 2 second window prior
% to stim presentation.
% Each row is a seperate individual stimulation presentation
%
% trajectoriesD - same as trajectoriesDDD but with just the raw diameter in
% microns
%
% trajectoryinfo - contains information about what happened in that
% particular stim; the rows correspond to the rows of trajectories.
% first column is the # of the movie
% second column is the vessel ID corresponding to that trajectory
% Third-fifth columns are the position of the vessel in the window (x,y,z
% in microns)
% sixth and seventh columns are the position of the stim on the screen
% (x,y as a fraction of the screen)
% eighth column is the global vessel ID which is a unique identifier for
% each vessel imaged in the whole window
% ninth column is the session#
%
% the above arrays are the results of averaging the stims in a given movie
% (i.e. if you present a particular stim 3 times in one movie, then you
% average over those 3 presentations). including a structure called
% "rawResults" which has the raw results of this without averaging,
% including the running data
%
% updates from v1.2
% reformat for processing the looped 2p aquisition (analog of
% poolStimData_v5_1.m)

function [trajectoriesDDD,trajectoryinfo,divingTracking,trajectoriesD,running,rawResults,trajectoriesGCaMP,pupilDisplacement,pupilRadius,stiminfotable] = poolStimData_diving_v1_4(regExpindicator,divingTracking,session,objectiveType,prestimTime)

trajectoriesGCaMP = [];         %this needs to be added in later
if nargin < 5 || isempty(prestimTime)
    prestimTime = 4.5;
end
fps = 10;
poststimTime = 9.5;
prestimBehavior = 3;
poststimBehavior = 3;
prestimBaseline = 2;
runningthresh = 300;     % exclude trajectories where mouse is running faster than this threshold
pupilthresh = 100;        % exclude trajectories where pupil is beyond this radius
minorAxisErrorThresh = .15;     %exclude points where the CI of the fit on the minor axis exceeds this threshold of the radius
fn = dir('*groupedDivingResults.mat');
doWithinMovieAverage = false;
if nargin < 4 || isempty(objectiveType)
    objectiveType = '25x';
end
if nargin < 3 || isempty(session)
    session = input('which imaging session is this: ');
end

if nargin < 2 || isempty(divingTracking)
    divingTracking = [];
    for n = 1:length(fn)
        load(fn(n).name,'results')
        if isempty(results(1).minoraxisFit)
            continue;
        end
        if contains(fn(n).name,regExpindicator{2})
            clear results
            continue;
        end
        for n2 = 1:length(results)
            numframes = size(results(n2).params,1);
            if size(results(n2).minoraxisFit,2) == numframes
                minoraxisFit = zeros(size(results(n2).minoraxisFit,2),size(results(n2).minoraxisFit,1),size(results(n2).minoraxisFit,3));
                for k1 = 1:size(minoraxisFit,2)
                    for k2 = 1:size(minoraxisFit,3)
                        minoraxisFit(:,k1,k2) = results(n2).minoraxisFit(k1,:,k2);
                    end
                end
                results(n2).minoraxisFit = minoraxisFit;
            end
            
            %get the corresponding stim metadata name:
            temp = strfind(fn(n).name,regExpindicator{1});
            temp = fn(n).name(temp:end);
            temp = regexp(temp,'\d*','Match');
            if isempty(temp)        % this means you didn't find any matching files
                continue;
            end
            strind = strfind(fn(n).name,regExpindicator{1});
            rootstring = [fn(n).name(1:(strind + length(regExpindicator{1}) - 1)),'_Metadata*'];
            vismetadata = dir([rootstring,temp{1},'.mat']);

            %% get the corresponding eye tracking 
            rootstringEye = fn(n).name(1:(strind + length(regExpindicator{1}) - 1));
            eyeTrackingFilename = dir([rootstringEye,temp{1},'*_EyeTracking.mat']);
            try
                [results(n2).eyeTracking,results(n2).eyeTrackingMovname] = getEyeTracking(eyeTrackingFilename(1).name);
            catch
                results(n2).eyeTracking = nan(size(results(n2).params,1)*size(results(n2).params,3),2);
                results(n2).eyeTrackingMovname = [];
                disp('cannot get eye tracking')
            end
            if size(results(n2).eyeTracking,1) < (size(results(n2).params,1)*size(results(n2).params,3))
                results(n2).eyeTracking = [nan((size(results(n2).params,1)*size(results(n2).params,3))-size(results(n2).eyeTracking,1),2);results(n2).eyeTracking];
            elseif size(results(n2).eyeTracking,1) > (size(results(n2).params,1)*size(results(n2).params,3))
                results(n2).eyeTracking = results(n2).eyeTracking(1:(size(results(n2).params,1)*size(results(n2).params,3)),:);
            end
            % reshape eyeTracking field to match the format of the other data
            results(n2).eyeTracking = permute(reshape(results(n2).eyeTracking,size(results(n2).params,1),size(results(n2).params,3),2),[2,1,3]);  
            [results(n2).stiminfo,results(n2).uniqueID,results(n2).stimMovieParams] = getStimInfo_v4_4(results(n2).moviename,vismetadata(1).name,results(n2));
            results(n2).running = results(n2).ballTracking;
            fieldsToRemove = {'stimID','ballTracking'};
            %results = rmfield(results,fieldsToRemove);
            results(n2).directory = cd;
            if exist('globalVesselID','var')
                results(n2).globalVesselID = globalVesselID;        %this value gives a unique identifier to each vessel imaged in the whole window; needs to be consistent across days. use the function "assignVesselID_manual_dev.m" currently to do so (11/20/2020)
            else
                results(n2).globalVesselID = [];
            end
            temp = results(n2);temp = rmfield(temp,fieldsToRemove);
            if ~isfield(temp,'objtype')
                temp.objtype = objectiveType;
            end
            divingTracking = [divingTracking;temp];
        end
        clear results validtracking globalVesslID
    end
end

if isempty(divingTracking)
    trajectoriesD = [];
    trajectoriesDDD = [];
    trajectoryinfo = [];
    running = [];
    return;
end

rawArray = [];      %this stores deltaD/D
rawArrayDiameter = [];  %this stores the raw diameter (in microns)
arrayInfo = [];
running = [];
pupilDisplacement = [];
pupilRadius = [];

for n = 1:length(divingTracking)       %iterate through each movie
    try
        meta = getSImetadata(divingTracking(n).moviename(1));
        curzoom = meta.zoom;
        if isfield(divingTracking(n),'objtype')
            micronperpixel = pixel2micron_2p(curzoom,divingTracking(n).objtype);
        else
            micronperpixel = pixel2micron_2p(curzoom,objectiveType);
        end
    catch
        curzoom = 3;
        disp('zoom error!')
    end
    
    %% figure out which one is the minor axis and get the error on that fit:
    numframes = max(size(divingTracking(n).running));
    if size(divingTracking(n).minoraxisFit,2) == numframes
        minoraxisFit  = zeros(size(divingTracking(n).minoraxisFit,2),size(divingTracking(n).minoraxisFit,1),size(divingTracking(n).minoraxisFit,3));
        for n2 = 1:size(minoraxisFit,3)
            minoraxisFit(:,:,n2) = divingTracking(n).minoraxisFit(:,:,n2)';
        end
        divingTracking(n).minoraxisFit = minoraxisFit;
    end
    rawTracking = divingTracking(n).minoraxisFit;
    for n2 = 1:size(divingTracking(n).params,3)
        if sum(isnan(divingTracking(n).minoraxisFit(:,1,n2)))/numel(isnan(divingTracking(n).minoraxisFit(:,1,n2))) > .25
            continue;
        end
        temp = divingTracking(n).params(:,:,n2);
        if size(divingTracking(n).ci,1) == 1        %if just doing the radon approach then no calculated ci's
            ci = zeros(size(divingTracking(n).params,1),1,size(divingTracking(n).params,3));
        else
            ci = divingTracking(n).ci(:,:,:,n2);
            ci = abs(diff(ci,1,3));
            temp(:,2) = temp(:,2).*temp(:,1);
            ci(:,2) = ci(:,2).*temp(:,2);           %again, update this if doing multiple ROIs in a single movie
            [~,I] = min(temp(:,1:2),[],2);
            I = sub2ind(size(ci),(1:size(ci,1))',I);
            ci = ci(I);
            ci = (ci')./temp;
        end
        
        %% not sure that this component is relevant anymore:
        temp = removeOutlierPoints(divingTracking(n).minoraxisFit(:,1,n2));
        temp(ci>minorAxisErrorThresh) = nan;
        rawTracking(:,1,n2) = temp;
    end
    
    vesselPositions = divingTracking(n).roiPosition;
    up = divingTracking(n).uniqueID;
    if isfield(divingTracking,'globalVesselID')
        globalVesselID = divingTracking(n).globalVesselID;
    else
        globalVesselID = [];
    end
    rn = divingTracking(n).running;
    
    numvessels = size(rawTracking,2);              %how was this working with this as the numvessels value??
    
    eyeTrack = divingTracking(n).eyeTracking;
    DDD = cell(1,numvessels);
    D = cell(1,numvessels);
    baselineFrames = ((divingTracking(n).preStimTime-prestimBaseline)*fps):(divingTracking(n).preStimTime*fps);
    for n2 = 1:numvessels
        ind = n2;
        [a,b] = averageForVessel(rawTracking,ind,baselineFrames);
        b = b*micronperpixel;
        DDD{n2} = permute(a,[3,1,2]);D{n2} = permute(b,[3,1,2]);
    end

    temp = divingTracking(n).stiminfo;    
    temp = temp(2,(divingTracking(n).preStimTime*fps+1),:);
    temp = permute(temp,[3,1,2]);
    curstimpos = zeros(size(temp,1),size(up,2));
    for n2 = 1:size(up,1)
        ind = find(temp == n2);
        curstimpos(ind,:) = repmat(up(n2,:),numel(ind),1);
    end
    
    for n2 = 1:size(DDD,1)
        if size(curstimpos,1) > size(DDD{n2},1)
            continue;
        end
        rawArray = [rawArray;DDD{n2}];
        rawArrayDiameter = [rawArrayDiameter;D{n2}];
        running = [running;rn];
        pupilRadius = [pupilRadius;eyeTrack(:,:,2)];
        pupilDisplacement = [pupilDisplacement;eyeTrack(:,:,1)];
        
        try
        temp = [repmat(n,size(D{n2},1),1),nan(size(D{n2},1),1),...
            repmat(vesselPositions(n2,:),size(D{n2},1),1),curstimpos];
        catch
            'a'
        end
        if isempty(globalVesselID)
            temp = [temp,nan(size(temp,1),1)];
        else
            temp = [temp,repmat(globalVesselID(uniqueVessels(n2)),size(temp,1),1)];
        end
        temp = [temp,repmat(session,size(temp,1),1)];
        arrayInfo = [arrayInfo;temp];
    end
end

trajectoriesDDD = rawArray;     %this stores deltaD/D trajectory
trajectoriesD = rawArrayDiameter;   %this stores just the diameter
trajectoryinfo = arrayInfo;
numstimboxes = round(size(arrayInfo,2)-7)/7;
% convert the array info into a table:
tableheader = {'movieNumber','localVesselID',...
    'vesselX_position','vesselY_position','vesselZ_position'};
tableheader = [tableheader,repmat({'stimStartX_position','stimStartY_position','stimMovieID',...
    'stimEndX_position','stimEndY_position','stimMovieSize','stimMovieDuration'},1,numstimboxes)];
tableheader = [tableheader,{'globalVesselID','sessionNumber'}];
stiminfotable = array2table(trajectoryinfo,'VariableNames',tableheader);

rawResults.trajectoriesDDD = trajectoriesDDD;
rawResults.trajectoriesD = trajectoriesD;
rawResults.trajectoryinfo = stiminfotable;
rawResults.running = running;
rawResults.pupilDisplacement = pupilDisplacement;
rawResults.pupilRadius = pupilRadius;

%% filter out trajectories that have too much running or pupil dilation:
stimduration = sum(divingTracking(1).stiminfo(3,:,1));
behaviorframes = ((divingTracking(1).preStimTime-prestimBehavior)*fps):((divingTracking(1).preStimTime+poststimBehavior)*fps+stimduration);
behaviorframes = behaviorframes(behaviorframes>0);
behaviorframes = behaviorframes(behaviorframes <= size(running,2));
avRun = nanmean(running(:,behaviorframes),2);
avPupil = nanmean(pupilRadius(:,behaviorframes),2);
i1 = find(avRun > runningthresh);
i2 = find(avPupil > pupilthresh);
ind = union(i1,i2);
trajectoriesDDD(ind,:) = nan;
trajectoriesD(ind,:) = nan;

%% average stims in given movies:
if doWithinMovieAverage
    temp = trajectoryinfo;
    [trajectoriesDDD,trajectoryinfo] = averageWithinMovie(trajectoriesDDD,trajectoryinfo);
    [trajectoriesD,~] = averageWithinMovie(trajectoriesD,trajectoryinfo);
    [running,~] = averageWithinMovie(running,temp);
    [pupilDisplacement,~] = averageWithinMovie(pupilDisplacement,temp);
    [pupilRadius,~] = averageWithinMovie(pupilRadius,temp);
    stiminfotable = array2table(trajectoryinfo,'VariableNames',tableheader);
end
%% crop to within the specificed prestim/poststim bounds:
validframes = (1+(divingTracking(1).preStimTime-prestimTime)*fps):((divingTracking(1).preStimTime+poststimTime)*fps+stimduration);
trajectoriesDDD = trajectoriesDDD(:,validframes);
trajectoriesD = trajectoriesD(:,validframes);
pupilRadius = pupilRadius(:,validframes);
pupilDisplacement = pupilDisplacement(:,validframes);
running = running(:,validframes);


%for a given trajectory, break it up into separate stims
% outpus: stimtraj is the broken out stimuli
% stimpos is an array of the positions of the stimulus on the monitor
% stimtrajRaw is just the raw tracking (not corrected for baseline
% diameter)
function [stimtraj,stimpos,stimtrajRaw,running] = breakupTraj(stiminfo,uniquePositions,traj,prestimFrames,poststimFrames,runningdata)

temp = [0,diff(stiminfo(3,:))];
stimstarts = find(temp == 1);
stimends = find(temp == -1)-1;

baseline = [];
stimtraj = [];
stimpos = [];
running = [];
if nargin > 5 && ~isempty(runningdata) && size(runningdata,2) == 1
    runningdata = runningdata';
end
    
for n = 1:length(stimstarts)
    curtraj = traj((stimstarts(n)-prestimFrames):(stimends(n)+poststimFrames));
    if nargin > 5 && ~isempty(runningdata)
        running = [running;runningdata((stimstarts(n)-prestimFrames):(stimends(n)+poststimFrames))];
    end
    baseline = [baseline;nanmean(curtraj(1:prestimFrames))];
    stimpos = [stimpos;uniquePositions(stiminfo(2,stimstarts(n)+1),:)];
    stimtraj = [stimtraj;curtraj];
end
stimtrajRaw = stimtraj;
stimtraj = stimtraj/mean(baseline);
return;

% since give multiple stims per movie, average all the stims within a movie
% with identical stim parameters
function [avTraj,avinfo] = averageIndividualMovies(traj,trajinfo)
indicator = trajinfo(:,[1,2,6,7]);
uniqueIndicator = unique(indicator,'rows');
avTraj = zeros(size(uniqueIndicator,1),size(traj,2));
avinfo = nan(size(uniqueIndicator,1),size(trajinfo,2));

for n = 1:size(uniqueIndicator,1)
    ind = ismember(indicator,uniqueIndicator(n,:),'rows');
    avTraj(n,:) = nanmean(traj(ind,:),1);
    avinfo(n,:) = mean(trajinfo(ind,:),1);
end

return;

% average all vessels within a given distance assuming identical stim
% parameters
function [avTraj,avinfo] = averagePositions(traj,trajinfo,distance)
indicator = trajinfo(:,3:5);

return;

% average all vessels with a stim position on the screen within a given
% distance (again keeping all other parameters)
function [avTraj,avinfo] = averageStimPositions(traj,info,distance)

return;


function r = getball(metadata,numframes)

load(metadata,'ballTrackingData')
if ~exist('ballTrackingData')
    r = [];
    disp('no running data available')
    return;
end
r = ballTrackingData{1};
r = diff(r(1:(length(r)-1)));
if numel(r) == numframes
    return;
else
    t = linspace(1,numframes,numel(r));
    r = interp1(1:numel(r),r,t);
end
return;


function [r,movname] = getEyeTracking(filename)
load(filename)
nf = size(center,1);
vi = find(~isnan(center(:,1)));
blanki = find(isnan(center(:,1)));
if numel(blanki)/nf > .7       %don't bother including movies where more than this threshold is non-valid
    r = nan(nf,2);
end
medpos = nanmedian(center,1);
center(:,1) = center(:,1)-medpos(1);
center(:,2) = center(:,2)-medpos(2);
center = sqrt(center(:,1).^2 + center(:,2).^2);
radius = radius(:,1);
center(blanki) = interp1(vi,center(vi),blanki);
radius(blanki) = interp1(vi,radius(vi),blanki);
blanki = find(radius<3); 
radius(blanki) = nan;
center(blanki) = nan;
r = [center,radius];
return;


function [stimtraj,stimtrajRaw] = averageForVessel(vesselTracking,ind,baselineFrames)
% this is taken from poolStimData_v5_1 for pial vessels where you have
% multiple kymographs per vessel. This is just trying to keep the same
% format but effectively all it's computing is deltaD/D
% given vessel. Outputs are stimtraj, which is deltaD/D and stimtrajRaw
% which is just diameter (in pixels)
baselineD = nanmean(vesselTracking(baselineFrames,ind,:),1);
stimtrajRaw = vesselTracking(:,ind,:);
repmat(baselineD,size(vesselTracking,1),1,1);
stimtraj = stimtrajRaw./repmat(baselineD,size(vesselTracking,1),1,1);

return;

function [trajAVG,trajinfoAVG] = averageWithinMovie(traj,trajinfo)
% averages identical stimuli within a given movie
sortinfo = trajinfo(:,[1;2;6;7;8;10]);
[~,ia,ib] = unique(sortinfo,'rows','stable');
trajinfoAVG = zeros(size(ia,1),size(trajinfo,2));
trajAVG = zeros(size(ia,1),size(traj,2));
for n = 1:size(ia,1)
    ind = find(ib == n);
    trajinfoAVG(n,:) = trajinfo(ia(n),:);
    trajAVG(n,:) = nanmean(traj(ind,:),1);
end
return;
